Sesión 6: Introducción a modelación
Jorge de la Vega
2023-02-23
En la sesión de hoy veremos:
- Uso de los modelos de regresión, caso particular de los {}. Los
modelos lineales en general abarcan:
- modelos de regresión (múltiple)
- modelos de diseños experimentales (variables y respuesta discretas)
- modelos para el análisis de covarianza (variables continuas y discretas, respuesta continua)
- los modelos lineales generalizados (regresión logística, binomial, Poisson)
- Trataremos de incluir:
- Cómo pensar sobre la formación de modelos
- media y varianza como funciones.
- estimación de parámetros
- diagnósticos
- inferencia
Regresión lineal simple
Los siguientes son datos del libro: “Understandable Statistics”. Los datos son para seguros de autos. Cada renglón es una región geográfica de Suecia.
library(readxl)
datos <- read_xls("../R_docentes/datos/slr06.xls")
head(datos)## # A tibble: 6 × 2
## X Y
## <dbl> <dbl>
## 1 108 392.
## 2 19 46.2
## 3 13 15.7
## 4 124 422.
## 5 40 119.
## 6 57 171.dim(datos) # dimensión del dataframe## [1] 63 2colnames(datos) <- c("n_claims","total_payment")Estadísticas descriptivas
summary(datos)## n_claims total_payment
## Min. : 0.0 Min. : 0.00
## 1st Qu.: 7.5 1st Qu.: 38.85
## Median : 14.0 Median : 73.40
## Mean : 22.9 Mean : 98.19
## 3rd Qu.: 29.0 3rd Qu.:140.00
## Max. :124.0 Max. :422.20cor(datos)## n_claims total_payment
## n_claims 1.0000000 0.9128782
## total_payment 0.9128782 1.0000000# Podemos generar una gráfica de la matriz de correlaciones con el paquete corrplot
library(corrplot)## corrplot 0.92 loadedcorrplot(cor(datos), method = "ellipse")
# gráfica de dispersión de puntos
plot(datos$n_claims, datos$total_payment)
abline(lm(total_payment ~ n_claims, data = datos)) # Agrega la línea de mínimos cuadrados
# regresión no paramétrica
a <- loess(total_payment ~ n_claims, data = datos) # Regresión no paramétrica
j <- order(a$x)
lines(a$x[j], a$fitted[j], col = "red", lwd=3)
# Versión ggplot:
library(ggplot2)
ggplot(datos,aes(n_claims, total_payment)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", se = F) +
geom_smooth(method = "loess", se = F, color = "red")## `geom_smooth()` using formula = 'y ~ x'## `geom_smooth()` using formula = 'y ~ x'
Ajuste del modelo de regresión
Ya lo hicimos arriba, para la gráfica, pero vamos a ver los detalles.
mod1 <- lm(total_payment ~ n_claims, data = datos)
mod1 # devuelve los parámetros estimados##
## Call:
## lm(formula = total_payment ~ n_claims, data = datos)
##
## Coefficients:
## (Intercept) n_claims
## 19.994 3.414summary(mod1) # nos da todos los detalles: estimación, errores estándar, anova. ##
## Call:
## lm(formula = total_payment ~ n_claims, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -86.561 -24.051 -0.347 23.432 83.977
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.9945 6.3678 3.14 0.0026 **
## n_claims 3.4138 0.1955 17.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 35.94 on 61 degrees of freedom
## Multiple R-squared: 0.8333, Adjusted R-squared: 0.8306
## F-statistic: 305 on 1 and 61 DF, p-value: < 2.2e-16model.matrix(mod1) # muestra la matriz de diseño X correspondiente al modelo## (Intercept) n_claims
## 1 1 108
## 2 1 19
## 3 1 13
## 4 1 124
## 5 1 40
## 6 1 57
## 7 1 23
## 8 1 14
## 9 1 45
## 10 1 10
## 11 1 5
## 12 1 48
## 13 1 11
## 14 1 23
## 15 1 7
## 16 1 2
## 17 1 24
## 18 1 6
## 19 1 3
## 20 1 23
## 21 1 6
## 22 1 9
## 23 1 9
## 24 1 3
## 25 1 29
## 26 1 7
## 27 1 4
## 28 1 20
## 29 1 7
## 30 1 4
## 31 1 0
## 32 1 25
## 33 1 6
## 34 1 5
## 35 1 22
## 36 1 11
## 37 1 61
## 38 1 12
## 39 1 4
## 40 1 16
## 41 1 13
## 42 1 60
## 43 1 41
## 44 1 37
## 45 1 55
## 46 1 41
## 47 1 11
## 48 1 27
## 49 1 8
## 50 1 3
## 51 1 17
## 52 1 13
## 53 1 13
## 54 1 15
## 55 1 8
## 56 1 29
## 57 1 30
## 58 1 24
## 59 1 9
## 60 1 31
## 61 1 14
## 62 1 53
## 63 1 26
## attr(,"assign")
## [1] 0 1Datos categóricos:
library(ggplot2)
datos2 <- iris[,c(1,5)] # Nos quedamos con las columnas uno y cinco
ggplot(datos2, aes(Sepal.Length)) +
geom_histogram(bins = 10) +
facet_wrap(vars(Species))
Estadísticas por categoría:
# tradicional
tapply(datos2$Sepal.Length,datos2$Species, summary)## $setosa
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.300 4.800 5.000 5.006 5.200 5.800
##
## $versicolor
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.900 5.600 5.900 5.936 6.300 7.000
##
## $virginica
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.900 6.225 6.500 6.588 6.900 7.900# tidy
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondatos2 %>%
group_by(Species) %>%
summarize(media = mean(Sepal.Length),
sd = sd(Sepal.Length),
n = n())## # A tibble: 3 × 4
## Species media sd n
## <fct> <dbl> <dbl> <int>
## 1 setosa 5.01 0.352 50
## 2 versicolor 5.94 0.516 50
## 3 virginica 6.59 0.636 50El modelo de regresión en este caso:
# Aquí R está considerando a Species como un factor, aunque es sólo un vector de caracteres.
# Se recomienda convertir Species a factor en caso de problemas.
mod2 <- lm(Sepal.Length ~ Species, data = datos2)
summary(mod2)##
## Call:
## lm(formula = Sepal.Length ~ Species, data = datos2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6880 -0.3285 -0.0060 0.3120 1.3120
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0060 0.0728 68.762 < 2e-16 ***
## Speciesversicolor 0.9300 0.1030 9.033 8.77e-16 ***
## Speciesvirginica 1.5820 0.1030 15.366 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5148 on 147 degrees of freedom
## Multiple R-squared: 0.6187, Adjusted R-squared: 0.6135
## F-statistic: 119.3 on 2 and 147 DF, p-value: < 2.2e-16X <- model.matrix(mod2)
head(X)## (Intercept) Speciesversicolor Speciesvirginica
## 1 1 0 0
## 2 1 0 0
## 3 1 0 0
## 4 1 0 0
## 5 1 0 0
## 6 1 0 0# También se puede quitar la ordenada del modelo
mod2a <- lm(Sepal.Length ~ Species + 0, data = datos2)
mod2a##
## Call:
## lm(formula = Sepal.Length ~ Species + 0, data = datos2)
##
## Coefficients:
## Speciessetosa Speciesversicolor Speciesvirginica
## 5.006 5.936 6.588Para extraer información del objeto
coefficients(mod1) # extrae los coeficientes## (Intercept) n_claims
## 19.994486 3.413824fitted(mod1) # extrae los valores ajustados## 1 2 3 4 5 6 7 8
## 388.68743 84.85713 64.37419 443.30861 156.54743 214.58243 98.51243 67.78802
## 9 10 11 12 13 14 15 16
## 173.61655 54.13272 37.06360 183.85802 57.54654 98.51243 43.89125 26.82213
## 17 18 19 20 21 22 23 24
## 101.92625 40.47743 30.23596 98.51243 40.47743 50.71890 50.71890 30.23596
## 25 26 27 28 29 30 31 32
## 118.99537 43.89125 33.64978 88.27096 43.89125 33.64978 19.99449 105.34007
## 33 34 35 36 37 38 39 40
## 40.47743 37.06360 95.09860 57.54654 228.23772 60.96037 33.64978 74.61566
## 41 42 43 44 45 46 47 48
## 64.37419 224.82390 159.96125 146.30596 207.75478 159.96125 57.54654 112.16772
## 49 50 51 52 53 54 55 56
## 47.30507 30.23596 78.02949 64.37419 64.37419 71.20184 47.30507 118.99537
## 57 58 59 60 61 62 63
## 122.40919 101.92625 50.71890 125.82302 67.78802 200.92713 108.75390residuals(mod1) # extrae los residuales del modelo (varios tipos de residual, con opción type = )## 1 2 3 4 5 6
## 3.8125698 -38.6571334 -48.6741920 -21.1086072 -37.1474282 -43.6824287
## 7 8 9 10 11 12
## -41.6124276 9.7119844 40.3834540 11.1672786 -16.1636036 64.2419834
## 13 14 15 16 17 18
## -34.0465449 -58.9124276 4.9087493 -20.2221329 32.9737488 10.4225729
## 19 20 21 22 23 24
## -25.8359564 14.4875724 -25.6774271 -2.0188978 1.3811022 -17.0359564
## 25 26 27 28 29 30
## -15.0953690 33.6087493 -21.8497800 9.8290430 -15.9912507 4.4502200
## 31 32 33 34 35 36
## -19.9944858 -36.1400748 -25.8774271 3.2363964 66.4013959 -0.3465449
## 37 38 39 40 41 42
## -10.6377229 -2.8603685 -21.0497800 -15.0156627 25.5258080 -22.4238994
## 43 44 45 46 47 48
## 21.3387483 6.4940425 -44.9547816 -86.5612517 -36.2465449 -19.5677219
## 49 50 51 52 53 54
## 28.7949258 9.6640436 64.0705137 28.6258080 -32.4741920 -39.1018392
## 55 56 57 58 59 60
## 8.2949258 14.3046310 72.0908074 35.9737488 36.6811022 83.9769839
## 61 62 63
## 27.7119844 43.6728656 78.7461017cooks.distance(mod1) # distancias de Cook para diagnósticos## 1 2 3 4 5 6
## 2.183354e-03 9.758274e-03 1.788235e-02 1.180318e-01 1.376184e-02 4.115139e-02
## 7 8 9 10 11 12
## 1.098506e-02 6.900603e-04 2.034982e-02 1.047074e-03 2.699222e-03 5.912148e-02
## 13 14 15 16 17 18
## 9.375163e-03 2.201761e-02 2.283703e-04 4.832807e-03 6.913315e-03 1.074464e-03
## 19 20 21 22 23 24
## 7.539119e-03 1.331520e-03 6.521462e-03 3.558483e-05 1.665282e-05 3.277967e-03
## 25 26 27 28 29 30
## 1.549092e-03 1.070539e-02 5.155686e-03 6.228258e-04 2.423610e-03 2.138727e-04
## 31 32 33 34 35 36
## 5.177366e-03 8.355639e-03 6.623448e-03 1.082146e-04 2.801473e-02 9.712964e-07
## 37 38 39 40 41 42
## 2.907237e-03 6.385946e-05 4.785060e-03 1.561876e-03 4.917973e-03 1.237114e-02
## 43 44 45 46 47 48
## 4.743961e-03 3.709949e-04 3.985862e-02 7.806390e-02 1.062591e-02 2.507445e-03
## 49 50 51 52 53 54
## 7.537632e-03 1.054847e-03 2.779156e-02 6.185042e-03 7.959831e-03 1.086948e-02
## 55 56 57 58 59 60
## 6.255010e-04 1.391051e-03 3.617134e-02 8.228508e-03 1.174684e-02 5.039858e-02
## 61 62 63
## 5.618318e-03 3.436581e-02 4.006284e-02# Un poco más avanzado:
library(broom)
tidy(mod1) # regresa las estadísticas en forma de dataframe para programación## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 20.0 6.37 3.14 2.60e- 3
## 2 n_claims 3.41 0.195 17.5 2.05e-25augment(mod1) # un dataframe con varias estadísticas asociadas a cada par de datos. ## # A tibble: 63 × 8
## total_payment n_claims .fitted .resid .hat .sigma .cooksd .std.resid
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 392. 108 389. 3.81 0.230 36.2 0.00218 0.121
## 2 46.2 19 84.9 -38.7 0.0163 35.9 0.00976 -1.08
## 3 15.7 13 64.4 -48.7 0.0188 35.7 0.0179 -1.37
## 4 422. 124 443. -21.1 0.318 36.1 0.118 -0.711
## 5 119. 40 157. -37.1 0.0245 35.9 0.0138 -1.05
## 6 171. 57 215. -43.7 0.0503 35.8 0.0412 -1.25
## 7 56.9 23 98.5 -41.6 0.0159 35.8 0.0110 -1.17
## 8 77.5 14 67.8 9.71 0.0182 36.2 0.000690 0.273
## 9 214 45 174. 40.4 0.0303 35.9 0.0203 1.14
## 10 65.3 10 54.1 11.2 0.0208 36.2 0.00105 0.314
## # … with 53 more rowsglance(mod1) # estadísticas generales del modelo## # A tibble: 1 × 12
## r.squared adj.r.squa…¹ sigma stati…² p.value df logLik AIC BIC devia…³
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.833 0.831 35.9 305. 2.05e-25 1 -314. 634. 640. 78797.
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## # variable names ¹adj.r.squared, ²statistic, ³deviancePredicciones
Usando el primer modelo con variables continuas. Necesitamos crear un arreglo con los valores de los predictores que queremos considerar en la evaluación
ncasos <- data.frame(n_claims=80:100)
predict(mod1,ncasos)## 1 2 3 4 5 6 7 8
## 293.1004 296.5142 299.9280 303.3418 306.7557 310.1695 313.5833 316.9971
## 9 10 11 12 13 14 15 16
## 320.4110 323.8248 327.2386 330.6524 334.0663 337.4801 340.8939 344.3077
## 17 18 19 20 21
## 347.7215 351.1354 354.5492 357.9630 361.3768# Serìa mejor tener en un sólo arreglo los valores nuevos y sus predicciones
predicciones <- ncasos %>%
mutate(total_payment = predict(mod1,ncasos))
# Podemos agregar las predicciones a una gráfica
plot(datos$n_claims, datos$total_payment)
abline(mod1)
points(predicciones, col = "red", pch = 19)
# Con ggplot
datos %>%
ggplot(aes(x = n_claims, y = total_payment)) +
geom_point() +
geom_smooth(method = "lm", se = F) +
geom_point(
data = predicciones,
color = "red"
)## `geom_smooth()` using formula = 'y ~ x'
También podemos obtener los intervalos de confianza para las predicciones. Hay dos opciones: una para la estimación de un valor promedio, y la otra para una nueva observación.
A <- as.data.frame(predict(mod1,interval = "confidence",level=0.8))
B <- as.data.frame(predict(mod1,interval = "prediction",level=0.8))## Warning in predict.lm(mod1, interval = "prediction", level = 0.8): predictions on current data refer to _future_ responses# añade intervalos de confianza a la gráfica
plot(datos$n_claims,datos$total_payment)
abline(mod1)
lines(datos$n_claims, A$lwr, col = "red")
lines(datos$n_claims, A$upr, col = "red")
lines(datos$n_claims, B$lwr, col = "green")
lines(datos$n_claims, B$upr, col = "green")
# Con ggplot:
datos %>%
ggplot(aes(x = n_claims, y = total_payment)) +
geom_point() +
geom_smooth(method = "lm", se = T) +
geom_line(aes(y=B$lwr), color = "orange", linetype = "dashed") +
geom_line(aes(y=B$upr), color = "orange", linetype = "dashed") ## `geom_smooth()` using formula = 'y ~ x'
Diagnósticos para el modelo
Hay una serie de gráficas para verificar los supuestos del modelo
par(mfrow = c(2,2)) # crea un arreglo para múltiples gráficas
plot(mod1) # gráficas de diagnóstico
Transformación de variables
Consideremos los siguientes datos: es un registro de 7 especies comunes de peces en ventas de mercado. Los datos son de Kaggle. Seleccionamos un grupo particular, el grupo de lubinas (perch)
Fish <- read.csv("https://raw.githubusercontent.com/jvega68/EA3/master/datos/Fish.csv")
unique(Fish$Species)## [1] "Bream" "Roach" "Whitefish" "Parkki" "Perch" "Pike"
## [7] "Smelt"perch <- Fish %>%
filter(Species == "Perch")
with(perch, plot(Length1,Weight)) # la relación no es lineal 
with(perch, plot(Length1^3,Weight)) # Vemos que la relación se lineariza con esta transformación
Para hacer la transformación en la regresión, podemos aplicar
funciones a los predictores, pero en el caso de
potencias, requerimos añadir el operador I(), porque el
símbolo de potencia ^ tiene una función especial, que se
verá más adelante
(mod3 <- lm(Weight ~ I(Length1^3), data = perch))##
## Call:
## lm(formula = Weight ~ I(Length1^3), data = perch)
##
## Coefficients:
## (Intercept) I(Length1^3)
## -0.1175 0.0168Para realizar la predicción, no es necesario elevar los datos al cubo, eso se hace como parte del modelo:
nuevos_datos <- data.frame( Length1 = seq(10,40,5))
datos_predicc <- data.frame(x = nuevos_datos,
y = predict(mod3, nuevos_datos))
with(perch, plot(Length1^3,Weight)) # Vemos que la relación se lineariza con esta transformación
abline(mod3)
# agrega los puntos rehaciendo la transformación para que los grafique correctamente.
points(datos_predicc$Length1^3, datos_predicc$y, col = "blue", pch = 19, cex = 2)
El tratamiento para transformar la respuesta es un poco diferente. Los siguientes datos también son de Kaggle y tienen información de tres campañas de marketing. Cada renglón es un anuncio.
Nos fijamos en tres variables: - spent: gasto en USD -
impressions: las veces que las personas vieron los anuncios
- clicks: las veces que las personas dieron clcks en los
anuncios
fb_ad <- read.csv("https://raw.githubusercontent.com/jvega68/EA3/master/datos/fb_adv_data.csv")
str(fb_ad)## 'data.frame': 1143 obs. of 15 variables:
## $ ad_id : int 708746 708749 708771 708815 708818 708820 708889 708895 708953 708958 ...
## $ reporting_start : chr "17/08/2017" "17/08/2017" "17/08/2017" "30/08/2017" ...
## $ reporting_end : chr "17/08/2017" "17/08/2017" "17/08/2017" "30/08/2017" ...
## $ campaign_id : chr "916" "916" "916" "916" ...
## $ fb_campaign_id : chr "103916" "103917" "103920" "103928" ...
## $ age : chr "30-34" "30-34" "30-34" "30-34" ...
## $ gender : chr "M" "M" "M" "M" ...
## $ interest1 : int 15 16 20 28 28 29 15 16 27 28 ...
## $ interest2 : int 17 19 25 32 33 30 16 20 31 32 ...
## $ interest3 : int 17 21 22 32 32 30 17 18 31 31 ...
## $ impressions : num 7350 17861 693 4259 4133 ...
## $ clicks : int 1 2 0 1 1 0 3 1 1 3 ...
## $ spent : num 1.43 1.82 0 1.25 1.29 ...
## $ total_conversion : int 2 2 1 1 1 1 1 1 1 1 ...
## $ approved_conversion: int 1 0 0 0 1 1 0 1 0 0 ...with(fb_ad, plot(spent, impressions)) # los puntos están muy concentrados en la parte baja de la gráfica
with(fb_ad, plot(sqrt(spent), sqrt(impressions))) # se aprecian mejor los datos de la parte baja
# filtramos los casos con 0 impresiones o con 0 gasto
fb1 <- fb_ad %>%
filter(!((spent < 10) | (impressions == 0)))
dim(fb1)## [1] 263 15with(fb1, plot(sqrt(spent), sqrt(impressions))) # se aprecian mejor los datos de la parte baja
mod4 <- lm(sqrt(impressions) ~ sqrt(spent), data = fb1)
mod4##
## Call:
## lm(formula = sqrt(impressions) ~ sqrt(spent), data = fb1)
##
## Coefficients:
## (Intercept) sqrt(spent)
## -32.06 65.98# prediccion
nvos_datos <- data.frame(spent = seq(0, 600, 100))
# la predicción devuelve los valores de la variable predictora en raíz cuadrada, por lo que necesitamos regresar
# a la escala original
datos_pred <- data.frame(
spent = seq(0, 600, 100),
sqrt_impressions = predict(mod4, nvos_datos))
datos_pred$impressions = datos_pred$sqrt_impressions^2 # agrega los datos en la escala original
plot(fb1$spent, fb1$impressions)
points(datos_pred$spent, datos_pred$impressions, col = "orange",pch = 19, cex = 2)
Regresión lineal múltiple
Revisaremos algunos puntos importantes en el análisis de regresión lineal múltiple, particularmente las diferencias de RLM con RLS. Muchas de las cosas que revisamos en RLS se aplican directamente en RLM.
Especificación de modelos en R en RLM
A continuación se muestra notación de modelos multivariados.
y ~ x\(y=b_0+b_1x\) Modelo de regresión lineal simpley ~ -1 + x\(y= b_1x\) Modelo sin ordenada al origeny ~ x + I(x^2)\(y=b_0 + b_1x+b_2x^2\) Modelo polinomial.y ~ x1 + x2\(y = b_0 + b_1x_1 + b_2x_2\) Modelo de primer orden sin interacción.y ~ x1:x2\(y = b_0 + b_1x_1x_2\) Modelo de interacción de primer ordeny ~ x1*x2\(y = b_0 + b_1x_1 + b_2x_2 + b_3x_1x_2\) Modelo de primer orden completo. Un modelo equivalente es: -y ~ x1 + x2 + x1:x2y ~ (x1 + x2 + x3)^2\(y = b_0 + b_1x_1 + b_2x_2 + b_3x_3 +b_4x_1x_2 +b_5x_2x_3 +b_6x_1x_3\) Interacciones dobles
Consideremos el siguiente ejemplo con el archivo de
Fish
Bream <- subset(Fish, Species == "Bream", select = -Species) # otra forma de tomar subconjuntos de datos.
mod5 <- lm(Weight ~ Length3 + Height + Width, data = Bream)
summary(mod5)##
## Call:
## lm(formula = Weight ~ Length3 + Height + Width, data = Bream)
##
## Residuals:
## Min 1Q Median 3Q Max
## -169.057 -26.519 -1.592 29.902 109.631
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1046.601 90.341 -11.585 8.54e-13 ***
## Length3 12.033 7.615 1.580 0.124248
## Height 62.820 16.685 3.765 0.000699 ***
## Width 45.898 35.415 1.296 0.204548
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 52.76 on 31 degrees of freedom
## Multiple R-squared: 0.942, Adjusted R-squared: 0.9364
## F-statistic: 167.8 on 3 and 31 DF, p-value: < 2.2e-16# Podemos aplicar directamente a los datos originales
lm(Weight ~ Length3 + Height + Width, data = Fish, subset = Species == "Bream")##
## Call:
## lm(formula = Weight ~ Length3 + Height + Width, data = Fish,
## subset = Species == "Bream")
##
## Coefficients:
## (Intercept) Length3 Height Width
## -1046.60 12.03 62.82 45.90lm(Weight ~ Length3 + Height + Width, data = Fish, subset = Species == "Roach")##
## Call:
## lm(formula = Weight ~ Length3 + Height + Width, data = Fish,
## subset = Species == "Roach")
##
## Coefficients:
## (Intercept) Length3 Height Width
## -319.924 8.608 -1.875 73.701Scatterplot matrix: datos en pares
plot(Fish[,-1], col = factor(Fish$Species))
# Equivalente en ggplot
library(GGally) # Perdón, omití decirles que instalaran este paquete## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2Fish %>%
group_by(Species) %>%
ggpairs(aes(color = Species), alpha = 0.4)## Warning in warn_if_args_exist(list(...)): Extra arguments: 'alpha' are being
## ignored. If these are meant to be aesthetics, submit them using the 'mapping'
## variable within ggpairs with ggplot2::aes or ggplot2::aes_string.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning in cor(x, y): the standard deviation is zero## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Otro paquete recién descubierto
library(scatterPlotMatrix)
Fish$Species <- factor(Fish$Species)
scatterPlotMatrix(Fish,zAxisDim = "Species")Gráfica de correlaciones:
corrplot(cor(Fish[,-1]), method = "ellipse", order = "hclust")
Para realizar selección de modelos de manera automática
library(MASS) # No requiere instalación, es de los paquetes básicos que se instalan junto con R##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectmod5 <- lm(Weight ~ (Length3 + Height + Width)^3, data = Bream)
mod5##
## Call:
## lm(formula = Weight ~ (Length3 + Height + Width)^3, data = Bream)
##
## Coefficients:
## (Intercept) Length3 Height
## 2593.952 -72.945 -384.968
## Width Length3:Height Length3:Width
## -121.926 11.110 2.541
## Height:Width Length3:Height:Width
## 47.000 -1.122mod6 <- stepAIC(mod5, # modelo inicial
scope = list(lower = ~1), # se puede dar un modelo upper y un lower
direction = "both" ) # foward o backward## Start: AIC=287.34
## Weight ~ (Length3 + Height + Width)^3
##
## Df Sum of Sq RSS AIC
## - Length3:Height:Width 1 1517.9 82985 285.99
## <none> 81467 287.34
##
## Step: AIC=285.99
## Weight ~ Length3 + Height + Width + Length3:Height + Length3:Width +
## Height:Width
##
## Df Sum of Sq RSS AIC
## - Height:Width 1 4.74 82990 283.99
## - Length3:Width 1 2103.13 85088 284.86
## - Length3:Height 1 2864.02 85849 285.18
## <none> 82985 285.99
##
## Step: AIC=283.99
## Weight ~ Length3 + Height + Width + Length3:Height + Length3:Width
##
## Df Sum of Sq RSS AIC
## - Length3:Width 1 2610.2 85600 283.07
## - Length3:Height 1 3145.5 86135 283.29
## <none> 82990 283.99
##
## Step: AIC=283.07
## Weight ~ Length3 + Height + Width + Length3:Height
##
## Df Sum of Sq RSS AIC
## - Length3:Height 1 705.7 86306 281.36
## - Width 1 4146.3 89746 282.73
## <none> 85600 283.07
##
## Step: AIC=281.36
## Weight ~ Length3 + Height + Width
##
## Df Sum of Sq RSS AIC
## - Width 1 4676 90982 281.21
## <none> 86306 281.36
## - Length3 1 6951 93256 282.07
## - Height 1 39467 125772 292.54
##
## Step: AIC=281.21
## Weight ~ Length3 + Height
##
## Df Sum of Sq RSS AIC
## <none> 90982 281.21
## - Length3 1 12717 103699 283.79
## - Height 1 62192 153173 297.44Comparación de modelos (siempre que estén anidados)
anova(mod5, mod6)## Analysis of Variance Table
##
## Model 1: Weight ~ (Length3 + Height + Width)^3
## Model 2: Weight ~ Length3 + Height
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 81467
## 2 32 90982 -5 -9514.6 0.6307 0.6779Transformación de Box-Cox: Determina la transformación optima de la variable de respuesta a normalidad.
boxcox(mod5)
Referencias
Buenas referencias para modelos lineales hay muchas, dentro de ellas:
- John Fox. Applied Regression Analysis and Generalized Linear Models, 3rd. Ed. Sage, 2016.
- John Fox & Sandford Weisberg. An R Companion to Applied
Regression, 2nd. Ed. Sage, 2011. Este libro tiene su propio
paquete:
car(companion to applied regression) - W. N. Venables & B. D. Ripley. Modern Applied Statistics
with S, 4th. Ed. Springer, 2002. Este libro es la base del paquete
MASS